/*
 * ActivityCoefficient.cpp
 *
 *  Created on: 23 May 2011
 *      Author: allan
 */

#include "ActivityCoefficient.h"

// C++ includes
#include <cmath>
#include <iostream>
#include <map>
#include <string>
using namespace std;

// FluidReact includes
#include "NumericalUtilities.h"

//----------------------------------------------------------------------------------------------------//
// Auxiliary Data and Functions
//----------------------------------------------------------------------------------------------------//
std::map<string,double> effectiveElectrostaticRadius = 
{
	{"H+",	  3.08},
	{"Li+",   1.64},
	{"Na+",   1.91},
	{"K+",    2.27},
	{"Rb+",   2.41},
	{"Cs+",   2.61},
	{"NH4+",  2.31},
	{"Ag+",   2.20},
	{"Au+",   2.31},
	{"Cu+",   1.90},
	{"Mg++",  2.54},
	{"Sr++",  3.00},
	{"Ca++",  2.87},
	{"Ba++",  3.22},
	{"Pb++",  3.08},
	{"Zn++",  2.62},
	{"Cu++",  2.60},
	{"Cd++",  2.85},
	{"Hg++",  2.98},
	{"Fe++",  2.62},
	{"Mn++",  2.68},
	{"Fe+++", 3.46},
	{"Al+++", 3.33},
	{"Au+++", 3.72},
	{"La+++", 3.96},
	{"Gd+++", 3.79},
	{"In+++", 3.63},
	{"Ca+++", 3.44},
	{"F-",    1.33},
	{"Cl-",   1.81},
	{"Br-",   1.96},
	{"I-",    2.20},
	{"OH-",   1.40},
	{"HS-",   1.84},
	{"NO3-",  2.81},
	{"HCO3-", 2.10},
	{"HSO4-", 2.37},
	{"ClO4-", 3.59},
	{"ReO4-", 4.23},
	{"SO4--", 3.15},
	{"CO3--", 2.81}
};

//std::map<string,double> effectiveElectrostaticRadius = 
//{
//	{"H[+]",	3.08},
//	{"Li[+]",   1.64},
//	{"Na[+]",   1.91},
//	{"K[+]",    2.27},
//	{"Rb[+]",   2.41},
//	{"Cs[+]",   2.61},
//	{"NH4[+]",  2.31},
//	{"Ag[+]",   2.20},
//	{"Au[+]",   2.31},
//	{"Cu[+]",   1.90},
//	{"Mg[2+]",  2.54},
//	{"Sr[2+]",  3.00},
//	{"Ca[2+]",  2.87},
//	{"Ba[2+]",  3.22},
//	{"Pb[2+]",  3.08},
//	{"Zn[2+]",  2.62},
//	{"Cu[2+]",  2.60},
//	{"Cd[2+]",  2.85},
//	{"Hg[2+]",  2.98},
//	{"Fe[2+]",  2.62},
//	{"Mn[2+]",  2.68},
//	{"Fe[3+]",  3.46},
//	{"Al[3+]",  3.33},
//	{"Au[3+]",  3.72},
//	{"La[3+]",  3.96},
//	{"Gd[3+]",  3.79},
//	{"In[3+]",  3.63},
//	{"Ca[3+]",  3.44},
//	{"F[-]",    1.33},
//	{"Cl[-]",   1.81},
//	{"Br[-]",   1.96},
//	{"I[-]",    2.20},
//	{"OH[-]",   1.40},
//	{"HS[-]",   1.84},
//	{"NO3[-]",  2.81},
//	{"HCO3[-]", 2.10},
//	{"HSO4[-]", 2.37},
//	{"ClO4[-]", 3.59},
//	{"ReO4[-]", 4.23},
//	{"SO4[2-]", 3.15},
//	{"CO3[2-]", 2.81}
//};

const double ElectrostaticParameterEta = 1.66027E+05;


const unsigned NumPointsOfDebyeHuckelA               = 8;
const unsigned NumPointsOfDebyeHuckelB               = 8;
const unsigned NumPointsOfDebyeHuckelParameterNaCl   = 12;
const unsigned NumPointsOfDebyeHuckelParameterNapClm = 12;

const double TemperatureOfDebyeHuckelA[]               = {0, 25.0, 60.0, 100.0, 150.0, 200.0, 250.0, 300.0};
const double TemperatureOfDebyeHuckelB[]               = {0, 25.0, 60.0, 100.0, 150.0, 200.0, 250.0, 300.0};
const double TemperatureOfDebyeHuckelParameterNaCl[]   = {25.0, 50.0, 75.0, 100.0, 125.0, 150.0, 175.0, 200.0, 225.0, 250.0, 275.0, 300.0};
const double TemperatureOfDebyeHuckelParameterNapClm[] = {25.0, 50.0, 75.0, 100.0, 125.0, 150.0, 175.0, 200.0, 225.0, 250.0, 275.0, 300.0};

const double ValuesOfDebyeHuckelA[]               = {0.4939, 0.5114, 0.5465, 0.5995, 0.6855, 0.7994, 0.9593, 1.2180};
const double ValuesOfDebyeHuckelB[]               = {0.3253, 0.3288, 0.3346, 0.6842, 0.3421, 0.3639, 0.3766, 0.3925};
const double ValuesOfDebyeHuckelParameterNaCl[]   = {2.47, 2.15, 1.79, 1.39, 0.93, 0.41, -0.18, -0.85, -1.64, -2.57, -3.71, -5.21};
const double ValuesOfDebyeHuckelParameterNapClm[] = {-9.77, -5.59, -2.43, 0.23, 2.44, 4.51, 6.38, 8.11, 9.73, 11.29, 12.71, 14.15}; 

const vector<double> EffectiveElectrostaticRadiiOfIons(const vector<string>& namesOfIons)
{
	vector<double> effRadii(namesOfIons.size());
	
	for(unsigned iIon = 0; iIon < namesOfIons.size(); ++iIon) 
	{
		auto iterator = effectiveElectrostaticRadius.find(namesOfIons[iIon]);
		
		if(iterator == effectiveElectrostaticRadius.end()) // TODO Improve the way this code handle this exception.
		{
			cerr << "There is no ion " << namesOfIons[iIon] << " in the database." << endl;
			exit(0);
		}
		
		effRadii[iIon] = iterator->second;
	}
	
	return effRadii;
}

const double EffectiveElectrostaticRadiusOfIon(string namesOfIon)
{
	double effRadius;
	
	auto iterator = effectiveElectrostaticRadius.find(namesOfIon);
	
	if(iterator == effectiveElectrostaticRadius.end()) // TODO Improve the way this code handle this exception.
	{
		cerr << "There is no ion " << namesOfIon << " in the database." << endl;
		exit(0);
	}
	
	effRadius = iterator->second;

	return effRadius;
}

const double DebyeHuckelParameterA(double T, double P)
{
	return LagrangeInterpolation(T, TemperatureOfDebyeHuckelA, ValuesOfDebyeHuckelA, NumPointsOfDebyeHuckelA);
}

const double DebyeHuckelParameterB(double T, double P)
{
	return LagrangeInterpolation(T, TemperatureOfDebyeHuckelB, ValuesOfDebyeHuckelB, NumPointsOfDebyeHuckelB);
}

const double HKFParameterNaCl(double T, double P)
{
	double R = 1.9858775;

	return 1.0E-3 * 1.0/(2.302585 * R * (T + 273.15)) * LagrangeInterpolation(T, TemperatureOfDebyeHuckelParameterNaCl, ValuesOfDebyeHuckelParameterNaCl, NumPointsOfDebyeHuckelParameterNaCl);
}

const double HKFParameterNapClm(double T, double P)
{
	return 1.0E-2 * LagrangeInterpolation(T, TemperatureOfDebyeHuckelParameterNapClm, ValuesOfDebyeHuckelParameterNapClm, NumPointsOfDebyeHuckelParameterNapClm);
}

const double ActivityCoefficientIdealModel()
{
	return 1.0;
}

const double ActivityCoefficientNeutralModelSetschenow(double I, double setschenow)
{
	return pow(10.0, setschenow * I);
}

const double ActivityCoefficientIonicModelHKF(unsigned iIon, const ParamsHKF& p)
{
	// the square root of the ionic strength
	const double sqrtI = sqrt(p.I);

	// the electrical charge of the current ion
	const int z = p.zOfIons[iIon];
	
	// the square of the electrical charge of the ion
	const int zSquared = z * z;
	
	// the effective electrostatic radius of the ion
	const double effRadius = p.effRadii[iIon];
	
	// the absolute Born coefficient of the ion
	const double omegaAbs = ElectrostaticParameterEta * zSquared / effRadius;

	// the Debye-Huckel ion size parameter of the ion as computed by Reed, 1982 and also in TOUGHREACT
	const double a = (z < 0) ? 2.0 * (effRadius + 1.91 * abs(z))/(abs(z) + 1.0) : 
			                   2.0 * (effRadius + 1.81 * abs(z))/(abs(z) + 1.0);
	
	// the activity coefficient of the ion in mole-fraction scale
	double activityCoeff = -(p.A * zSquared * sqrtI)/(1.0 + a * p.B * sqrtI) + (omegaAbs * p.bNaCl + p.bNapClm - 0.19 * (abs(z) - 1.0)) * p.I;

	return pow(10.0, activityCoeff);
}

const double OsmoticCoefficientModelHKF(const ParamsHKF& p)
{
	// the number of ions
	const unsigned numOfIons = p.effRadii.size();
	
	// the square root of the ionic strength
	const double sqrtI = sqrt(p.I);

	// the Born coefficient of the ion H[+]
	const double omegaH = 0.5387E+05;
	
	// the osmotic coefficient of the aqueous phase
	double osmoticCoefficient = 0.0;
	
	// iterate over all ions in the aqueous phase	
	for(unsigned iIon = 0; iIon < numOfIons; ++iIon)
	{
		// the electrical charge of the current ion species
		const int z = p.zOfIons[iIon];
		
		// the square of the electrical charge
		const int zSquared = z * z;
		
		// the effective electrostatic radius of the current ion
		const double effRadius = p.effRadii[iIon];
		
		// the Born coefficient of the current ion
		const double omega = ElectrostaticParameterEta * zSquared / effRadius - z * omegaH;
		
		// the Debye-Huckel ion size parameter of the ion as computed by Reed, 1982 and also in TOUGHREACT
		const double a = (z < 0) ? 2.0 * (effRadius + 1.91 * abs(z))/(abs(z) + 1.0) : 
				                   2.0 * (effRadius + 1.81 * abs(z))/(abs(z) + 1.0);
		
		// the Lambda parameter of the current ion
		const double Lambda = 1.0 + a * p.B * sqrtI;
		
		// the sigma parameter of the current ion
		const double sigma = 3.0/pow(a * p.B * sqrtI, 3) * (Lambda - 1.0/Lambda - 2.0 * log(Lambda));
		
		// calculate the psi contribution of the current ion in mole-fraction scale
		double psi = p.A * zSquared * sqrtI * sigma/3.0 - 0.5 * (omega * p.bNaCl + p.bNapClm - 0.19 * (abs(z) - 1.0)) * p.I;
		
		osmoticCoefficient += p.mTi[iIon] * psi;
	}
	
	return -2.303/p.mT * osmoticCoefficient;
}

const double ActivityWaterModelHKF(const ParamsHKF& p)
{
	// the osmotic coefficient of the aqueous phase
	const double osmoticCoefficient = OsmoticCoefficientModelHKF(p);
	
	// the activity of water  in mole-fraction scale
	const double activityWater = exp(-osmoticCoefficient * p.mT / 55.508);
	
	return activityWater;
}

